## Deserving Government Assistance? 
## Public Support for Aid to Struggling Firms and Workers
#############################################################
## Chris Witko
## Tobias Heinrich
#############################################################
## 04) Estimate average Treatment Moderation Effects
#############################################################


## Load prepped data
####################
data <- read.csv(file="output/Prepped_survey_data.csv")
CCES <- read.csv(file="output/Prepped_CCES_data.csv")

translate_reason <- data.frame(reason=c("Laid off due to domestic competition", "Laid off due to foreign competition",
                                        "Operations were moved to a foreign country", 
                                        "Laid off due to company mismanagement", "Fired for poor performance", 
                                        "Job was replaced by technology (like robot, artificial intelligence)",
                                        "Laid off due to economic recession", "Company was closed due to pandemic public health restrictions",
                                        "Economic recession", "Competition due to superior use of technology (like robots, artificial intelligence)",
                                        "Mismanagement", "Foreign competition due to lower labor cost",
                                        "Operations closed due to pandemic public health restrictions"),
                               prettyreason=c("Domestic competition", "Foreign competition", "Offshoring",  
                                              "Company mismanagement", "Poor performance",
                                              "Technology", "Recession", "Public health",
                                              "Recession", "Technology", "Mismanagement", "Foreign competiton", "Public health"))
data$W_reasonforjobloss2 <- mapvalues(data$W_reasonforjobloss,
                                      from=translate_reason$reason,
                                      to=translate_reason$prettyreason)
data$F_reasonforstruggle2 <- mapvalues(data$F_reasonforstruggle,
                                       from=translate_reason$reason,
                                       to=translate_reason$prettyreason)

## Ebalance setup
CCES$id <- sample(1:10, size=nrow(CCES), replace=TRUE)
eb_form <- id ~ 1 + Gender_female + Age + Education_uni + I(Party == "Democrat") + I(Party == "Republican") +
  + I(Job == "Full-time") + I(Job == "Unemployed")
X_cces <- model.matrix(eb_form, CCES)[, -1]

## Formula for model; demographics only in ATME
form <- Choice ~ 1 + I(Inequality == "Very") + I(Inequality == "Somewhat") + Age + Gender_female + 
  + I(RaceEthnicity == "White") + I(RaceEthnicity == "Black") + I(RaceEthnicity == "Hispanic") +
  + I(Ideology == "Conservative") + I(Ideology == "Liberal") +
  + I(Ideology == "NS") + I(Party == "Democrat") + I(Party == "Republican") + 
  + I(Income == "Less than $30,000") + I(Income == "$60,000 - $89,999") + I(Income == "$90,000 - $119,999")  +
  + I(Income == "$120,000 or more") + I(Class == "Lower") + I(Class == "Upper") + I(Class == "Working") +
  + Education_uni + I(Job == "Full-time") + I(Job == "Student") +  I(Job == "Unemployed")

## Storage
all_ATMEs <- vector("list", 3)


## Looping through weighted/ unweighted models
for(j in 1:2)
{
  ## ATME for Hypothesis WM1
  ## Worker's salary
  ##########################
  ## High salary worker
  this_data <- subset(data, Which_experiment == "Worker" & W_salary == "$160,000")
  X_data <- model.matrix(eb_form, this_data)[, -1]
  X_eb <- rbind(X_data, X_cces)
  eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
  eb <- ebalance.trim(eb)
  w <- eb$w / sum(eb$w) * nrow(this_data)
  if(j == 1) mod_T1 <- lm(form, data=this_data, weights=w)
  if(j == 2) mod_T1 <- lm(form, data=this_data)
  VCV_T1 <- vcovBS(x=mod_T1, cluster=this_data$id)
  beta_T1 <- rmvnorm(1000, coef(mod_T1), VCV_T1)
  
  ## Low salary worker
  this_data <- subset(data, Which_experiment == "Worker" & W_salary == "$20,000")
  X_data <- model.matrix(eb_form, this_data)[, -1]
  X_eb <- rbind(X_data, X_cces)
  eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
  eb <- ebalance.trim(eb)
  w <- eb$w / sum(eb$w) * nrow(this_data)
  if(j == 1) mod_T0 <- lm(form, data=this_data, weights=w)
  if(j == 2) mod_T0 <- lm(form, data=this_data)
  VCV_T0 <- vcovBS(x=mod_T0, cluster=this_data$id)
  beta_T0 <- rmvnorm(1000, coef(mod_T0), VCV_T0)
  
  ## Package up
  all_ATMEs[[1]] <- list(what="ATME for Hypothesis WM1",
                         mod_T1=mod_T1, VCV_T1=VCV_T1, mod_T0=mod_T0, VCV_T0=VCV_T0,
                         ATME_ineq=beta_T1[, 2] - beta_T0[, 2],
                         ATME_party=(rowSums(beta_T1[, c(10, 12)]) - rowSums(beta_T1[, c(9, 13)])) - 
                           (rowSums(beta_T0[, c(10, 12)]) - rowSums(beta_T0[, c(9, 13)])))
  
  
  ## ATME for Hypothesis WF1
  ## Typical worker's salary
  ##########################
  ## High salary worker
  this_data <- subset(data, Which_experiment == "Firm" & F_typicalworkersalary == "$160,000")
  X_data <- model.matrix(eb_form, this_data)[, -1]
  X_eb <- rbind(X_data, X_cces)
  eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
  eb <- ebalance.trim(eb)
  w <- eb$w / sum(eb$w) * nrow(this_data)
  if(j == 1) mod_T1 <- lm(form, data=this_data, weights=w)
  if(j == 2) mod_T1 <- lm(form, data=this_data)
  VCV_T1 <- vcovBS(x=mod_T1, cluster=this_data$id)
  beta_T1 <- rmvnorm(1000, coef(mod_T1), VCV_T1)
  
  ## Low salary worker
  this_data <- subset(data, Which_experiment == "Firm" & F_typicalworkersalary == "$20,000")
  X_data <- model.matrix(eb_form, this_data)[, -1]
  X_eb <- rbind(X_data, X_cces)
  eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
  eb <- ebalance.trim(eb)
  w <- eb$w / sum(eb$w) * nrow(this_data)
  if(j == 1) mod_T0 <- lm(form, data=this_data, weights=w)
  if(j == 2) mod_T0 <- lm(form, data=this_data)
  VCV_T0 <- vcovBS(x=mod_T0, cluster=this_data$id)
  beta_T0 <- rmvnorm(1000, coef(mod_T0), VCV_T0)
  
  ## Package up
  all_ATMEs[[2]] <- list(what="ATME for Hypothesis WF1",
                         mod_T1=mod_T1, VCV_T1=VCV_T1, mod_T0=mod_T0, VCV_T0=VCV_T0,
                         ATME_ineq=beta_T1[, 2] - beta_T0[, 2],
                         ATME_party=(rowSums(beta_T1[, c(10, 12)]) - rowSums(beta_T1[, c(9, 13)])) - 
                           (rowSums(beta_T0[, c(10, 12)]) - rowSums(beta_T0[, c(9, 13)])))
  
  
  
  ## ATME for Hypothesis WF2
  ## CEO salary
  ##########################
  ## High salary CEO
  this_data <- subset(data, Which_experiment == "Firm" & F_ceopay == "$10 million")
  X_data <- model.matrix(eb_form, this_data)[, -1]
  X_eb <- rbind(X_data, X_cces)
  eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
  eb <- ebalance.trim(eb)
  w <- eb$w / sum(eb$w) * nrow(this_data)
  if(j == 1) mod_T1 <- lm(form, data=this_data, weights=w)
  if(j == 2) mod_T1 <- lm(form, data=this_data)
  VCV_T1 <- vcovBS(x=mod_T1, cluster=this_data$id)
  beta_T1 <- rmvnorm(1000, coef(mod_T1), VCV_T1)
  
  ## Low salary CEO
  this_data <- subset(data, Which_experiment == "Firm" & F_ceopay == "$200,000")
  X_data <- model.matrix(eb_form, this_data)[, -1]
  X_eb <- rbind(X_data, X_cces)
  eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
  eb <- ebalance.trim(eb)
  w <- eb$w / sum(eb$w) * nrow(this_data)
  if(j == 1) mod_T0 <- lm(form, data=this_data, weights=w)
  if(j == 2) mod_T0 <- lm(form, data=this_data)
  VCV_T0 <- vcovBS(x=mod_T0, cluster=this_data$id)
  beta_T0 <- rmvnorm(1000, coef(mod_T0), VCV_T0)
  
  ## Package up
  all_ATMEs[[3]] <- list(what="ATME for Hypothesis WF2",
                         mod_T1=mod_T1, VCV_T1=VCV_T1, mod_T0=mod_T0, VCV_T0=VCV_T0,
                         ATME_ineq=beta_T1[, 2] - beta_T0[, 2],
                         ATME_party=(rowSums(beta_T1[, c(10, 12)]) - rowSums(beta_T1[, c(9, 13)])) - 
                           (rowSums(beta_T0[, c(10, 12)]) - rowSums(beta_T0[, c(9, 13)])))
  
  ## Graph for ATMEs
  ##################
  ## Loop through inequality/ party moderation
  for(p in 1:2)
  {
    all_ATMEs2 <- c()
    for(i in 1:length(all_ATMEs))
    {
      if(p == 1)
      {
        all_ATMEs2 <- rbind(all_ATMEs2, 
                            data.frame(what=all_ATMEs[[i]]$what,
                                       Mean=mean(all_ATMEs[[i]]$ATME_ineq), Q025=quantile(all_ATMEs[[i]]$ATME_ineq, 0.025),
                                       Q975=quantile(all_ATMEs[[i]]$ATME_ineq, 0.975), PP=mean(all_ATMEs[[i]]$ATME_ineq > 0)))
      }
      if(p == 2)
      {
        all_ATMEs2 <- rbind(all_ATMEs2, 
                            data.frame(what=all_ATMEs[[i]]$what,
                                       Mean=mean(all_ATMEs[[i]]$ATME_party), Q025=quantile(all_ATMEs[[i]]$ATME_party, 0.025),
                                       Q975=quantile(all_ATMEs[[i]]$ATME_party, 0.975), PP=mean(all_ATMEs[[i]]$ATME_party > 0)))
      }
    }
    all_ATMEs2$Set <- ifelse(grepl(x=all_ATMEs2$what, pattern="WF") == TRUE, "Firms", "Workers")
    all_ATMEs2$Change <- NA
    all_ATMEs2$Change[grepl(x=all_ATMEs2$what, pattern="WM1")] <- "Worker salary\n$20k -> $160k"
    all_ATMEs2$Change[grepl(x=all_ATMEs2$what, pattern="WF1")] <- "Typical worker salary\n$20k -> $160k"
    all_ATMEs2$Change[grepl(x=all_ATMEs2$what, pattern="WF2")] <- "CEO pay\n$200k -> $10m"
    all_ATMEs2[, c("Mean", "Q025", "Q975")] <- 100 * all_ATMEs2[, c("Mean", "Q025", "Q975")]
    all_ATMEs2$Sig <- ifelse(all_ATMEs2$PP > 0.975 | all_ATMEs2$PP <= 0.025, "Sig", "Insig")
    all_ATMEs2$Mean_label <- format(round(all_ATMEs2$Mean, di=0), nsmall=0)
    
    g <- ggplot(data=all_ATMEs2, aes(x=Change, y=Mean))
    g <- g + geom_hline(yintercept=0, linewidth=0.8) + geom_point(size=0.002)
    g <- g + coord_flip() + theme_minimal()
    g <- g + geom_linerange(data=subset(all_ATMEs2, Sig == "Sig"), colour="black", 
                            aes(x=Change, ymin=Q025, ymax=Q975), linewidth=.99)
    g <- g + geom_linerange(data=subset(all_ATMEs2, Sig == "Insig"), colour="gray60", 
                            aes(x=Change, ymin=Q025, ymax=Q975), linewidth=.99)
    g <- g + geom_point(data=subset(all_ATMEs2, Sig == "Insig"), size=3.5, colour="gray60")
    g <- g + geom_label(data=subset(all_ATMEs2, Sig == "Sig"),
                        aes(label=Mean_label), size=3.6)
    g <- g + xlab("") + ylab(paste0("Average treatment moderation effect due to ", ifelse(p == 1, "inequality", "party affiliation")))
    g <- g + facet_wrap(~ Set, ncol=1, scales="free_y")
    g <- g + theme(axis.text.y=element_text(hjust=1, size=rel(1.34)),
                   axis.text.x=element_text(size=rel(1.34)),
                   strip.text=element_text(hjust=0, size=rel(1.2)),
                   plot.title=element_text(size=rel(1.44), hjust=0), 
                   legend.position="none") 
    if(p == 2 & j == 1) g <- g + ggtitle("Using reweighted models")
    if(p == 2 & j == 2) g <- g + ggtitle("Not using reweighted models")
    ggsave(file=paste0("output/ATMEs-", ifelse(p == 1, "inequality", "party"), "-", 
                       ifelse(j == 1, "weighted", "unweighted"), ".pdf"), 
           width=10, height=5.4 + ifelse(j == 1 & p == 1, 0, -1.2), device=cairo_pdf)
  }
  
  ## Big table for all constituent models
  #######################################
  all_coefs <- rbind(data.frame(What=all_ATMEs[[1]]$what, Subset="High salary",
                                Names=names(coef(all_ATMEs[[1]]$mod_T1)),
                                Coef=coef(all_ATMEs[[1]]$mod_T1), SE=sqrt(diag(all_ATMEs[[1]]$VCV_T1)),
                                N=length(all_ATMEs[[1]]$mod_T1$residuals),
                                Sigma=summary(all_ATMEs[[1]]$mod_T1)$sigma),
                     data.frame(What=all_ATMEs[[1]]$what, Subset="Low salary",
                                Names=names(coef(all_ATMEs[[1]]$mod_T0)),
                                Coef=coef(all_ATMEs[[1]]$mod_T0), SE=sqrt(diag(all_ATMEs[[1]]$VCV_T0)),
                                N=length(all_ATMEs[[1]]$mod_T0$residuals),
                                Sigma=summary(all_ATMEs[[1]]$mod_T0)$sigma),
                     data.frame(What=all_ATMEs[[2]]$what, Subset="High salary",
                                Names=names(coef(all_ATMEs[[2]]$mod_T1)),
                                Coef=coef(all_ATMEs[[2]]$mod_T1), SE=sqrt(diag(all_ATMEs[[2]]$VCV_T1)),
                                N=length(all_ATMEs[[2]]$mod_T1$residuals),
                                Sigma=summary(all_ATMEs[[2]]$mod_T1)$sigma),
                     data.frame(What=all_ATMEs[[2]]$what, Subset="Low salary",
                                Names=names(coef(all_ATMEs[[2]]$mod_T0)),
                                Coef=coef(all_ATMEs[[2]]$mod_T0), SE=sqrt(diag(all_ATMEs[[2]]$VCV_T0)),
                                N=length(all_ATMEs[[2]]$mod_T0$residuals),
                                Sigma=summary(all_ATMEs[[2]]$mod_T0)$sigma),
                     data.frame(What=all_ATMEs[[3]]$what, Subset="High CEO pay",
                                Names=names(coef(all_ATMEs[[3]]$mod_T1)),
                                Coef=coef(all_ATMEs[[3]]$mod_T1), SE=sqrt(diag(all_ATMEs[[3]]$VCV_T1)),
                                N=length(all_ATMEs[[3]]$mod_T1$residuals),
                                Sigma=summary(all_ATMEs[[3]]$mod_T1)$sigma),
                     data.frame(What=all_ATMEs[[3]]$what, Subset="Low CEO pay",
                                Names=names(coef(all_ATMEs[[3]]$mod_T0)),
                                Coef=coef(all_ATMEs[[3]]$mod_T0), SE=sqrt(diag(all_ATMEs[[3]]$VCV_T0)),
                                N=length(all_ATMEs[[3]]$mod_T0$residuals),
                                Sigma=summary(all_ATMEs[[3]]$mod_T0)$sigma))
  rownames(all_coefs) <- NULL
  all_coefs$Set <- ifelse(grepl(x=all_coefs$What, pattern="WF") == TRUE, "Firms", "Workers")
  all_coefs$order <- NA
  all_coefs$order[all_coefs$Subset == "Low CEO pay"] <- 5
  all_coefs$order[all_coefs$Subset == "High CEO pay"] <- 6
  all_coefs$order[all_coefs$Subset == "Low salary" & all_coefs$Set == "Workers"] <- 1
  all_coefs$order[all_coefs$Subset == "High salary" & all_coefs$Set == "Workers"] <- 2
  all_coefs$order[all_coefs$Subset == "Low salary" & all_coefs$Set == "Firms"] <- 3
  all_coefs$order[all_coefs$Subset == "High salary" & all_coefs$Set == "Firms"] <- 4
  
  
  ## Prettify names
  all_coefs$PrettyNames <- mapvalues(x=all_coefs$Names,
                                     from=c("(Intercept)", "I(Inequality == \"Very\")TRUE",
                                            "I(Inequality == \"Somewhat\")TRUE",
                                            "Age", "Gender_female", "I(RaceEthnicity == \"White\")TRUE",
                                            "I(RaceEthnicity == \"Black\")TRUE", "I(RaceEthnicity == \"Hispanic\")TRUE",
                                            "I(Ideology == \"Conservative\")TRUE", "I(Ideology == \"Liberal\")TRUE",
                                            "I(Ideology == \"NS\")TRUE", "I(Party == \"Democrat\")TRUE",
                                            "I(Party == \"Republican\")TRUE", "I(Income == \"Less than $30,000\")TRUE",
                                            "I(Income == \"$60,000 - $89,999\")TRUE", "I(Income == \"$90,000 - $119,999\")TRUE",
                                            "I(Income == \"$120,000 or more\")TRUE",
                                            "I(Class == \"Lower\")TRUE", "I(Class == \"Upper\")TRUE", "I(Class == \"Working\")TRUE",
                                            "Education_uni", "I(Job == \"Full-time\")TRUE",
                                            "I(Job == \"Student\")TRUE", "I(Job == \"Unemployed\")TRUE"),
                                     to=c("Intercept", "Inequality, very", "Inequality, somewhat",
                                          "Age", "Gender, female", "Race/ ethnicity, White", "Race/ ethnicity, Black",
                                          "Race/ ethnicity, Hispanic", "Ideology, conservative", "Ideology, liberal",
                                          "Ideology, not sure", "Party, Democrat", "Party, Republican",
                                          "Income, less than 30k", "Income, 60-89k", "Income, 90-119k",
                                          "Income, more than 120k", "Class, lower", "Class, upper", "Class, working",
                                          "Education, university", "Job, full-time", "Job, student", "Job, unemployed"))
  
  ## Make table
  u_vars <- c("Inequality, very", "Inequality, somewhat",
              "Age", "Gender, female", "Race/ ethnicity, White", "Race/ ethnicity, Black",
              "Race/ ethnicity, Hispanic", "Ideology, conservative", "Ideology, liberal",
              "Ideology, not sure", "Party, Democrat", "Party, Republican",
              "Income, less than 30k", "Income, 60-89k", "Income, 90-119k",
              "Income, more than 120k", "Class, lower", "Class, upper", "Class, working",
              "Education, university", "Job, full-time", "Job, student", "Job, unemployed", "Intercept")
  se_mat <- m_mat <- matrix("", length(u_vars), 6)
  T_vec <- N_vec <- sigma_vec <- rep(NA, 6)
  for(i in 1:6)
  {
    tmp <- subset(all_coefs, order == i)
    N_vec[i] <- tmp$N[1]
    sigma_vec[i] <- format(round(tmp$Sigma[1], di=2), nsmall=2)
    T_vec[i] <- tmp$Subset[1]
    for(k in 1:length(u_vars))
    {
      tmp2 <- subset(tmp, PrettyNames == u_vars[k])
      m_mat[k, i] <- format(round(tmp2$Coef, di=2), nsmall=2)
      se_mat[k, i] <- paste0("(", format(round(tmp2$SE, di=2), nsmall=2), ")")
    }
  }
  
  T_vec <- str_replace_all(string=T_vec, pattern="salary", replacement="pay")
  T_vec <- str_replace_all(string=T_vec, pattern="Low ", replacement="L ")
  T_vec <- str_replace_all(string=T_vec, pattern="High ", replacement="H ")
  s1 <- c(rep("Workers", 2), rep("Firms", 4))
  tab <- c("\\begin{table}[ht!]", "\\begin{center}", "\\scriptsize", 
           "\\renewcommand{\\arraystretch}{1.00}",
           paste0(c("\\begin{tabular}{@{\\extracolsep{4pt}}l", paste0(rep("c", 6)), "} "), collapse=""),     
           "\\\\[-1.8ex]\\hline \\\\[-1.8ex]")
  tab <- c(tab,
           paste0(c(" ", paste0("\\footnotesize{\\textbf{ ", s1, "}}")), collapse=" & "), " \\\\", 
           "[-1.8ex] \\\\ \\hline",
           paste0(c(" ", paste0("\\footnotesize{\\textbf{ ", T_vec, "}}")), collapse=" & "), " \\\\")
  
  for(k in 1:nrow(m_mat))
  {
    tab <- c(tab, 
             paste0(u_vars[k], " & ", paste0(m_mat[k, ], collapse=" & "), "\\\\"),
             paste0(" & ", paste0(se_mat[k, ], collapse=" & "), "\\\\"))
  }
  tab <- c(tab,
           paste0("Residual SE", " & ", paste0(sigma_vec, collapse=" & "), "\\\\"),
           paste0("Observations", " & ", paste0(N_vec, collapse=" & "), "\\\\"),
           "\\end{tabular}", 
           "\\caption{\\textbf{Coefficient estimates for ATME models.} 'L pay' and 'H pay' refer to low and high pay subsets of profile realizations for workers and firms; 'L CEO pay' and 'H CEO pay' indicate models estimated on subset with low and high CEO pay. Standard errors in parentheses.}",
           paste0("\\label{Tab:ATMEs", j, "}"),
           "\\end{center}", "\\end{table}")
  write(tab, file=paste0("output/Table-ATME-", ifelse(j == 1, "weighted", "unweighted"), ".tex"))
}


## Garbage collection
#####################
rm(list=ls())

